### Importing the dataset

fire_all <- read_csv("Data/Incidents_Responded_to_by_Fire_Companies.csv") 
fire_all$year <- substr(fire_all$INCIDENT_DATE_TIME, 7, 10)
fire <- fire_all%>% 
  filter(HIGHEST_LEVEL_DESC == "7 - Signal 7-5") %>%
  filter(year == 2015) %>%
  select(-Latitude, -Longitude)
save(fire, file = "Data/Fire.rda")
load("Data/Fire.rda")
### Merging in Latitude and Longitude
latlong <- read_csv("data/severe_incidents.csv") %>%
  select(IM_INCIDENT_KEY, Latitude, Longitude)

Merged <- left_join(fire, latlong, by = "IM_INCIDENT_KEY")
### Importing the FireHouse dataset

Firehouses <- read_csv("data/FDNY_Firehouse_Listing.csv") %>%
  select(FacilityName, Borough, Latitude, Longitude)

Question 1

content <- paste("Type of incident:",Merged$INCIDENT_TYPE_DESC,"<br/>",
                 "Action:",Merged$ACTION_TAKEN1_DESC,"<br/>",
                 "Detector Status:", Merged$DETECTOR_PRESENCE_DESC,"<br/>")




FireMap1 <- leaflet(Merged) %>%
  addTiles() %>%    # Add OpenStreetMap map tiles
  addCircles(lng = ~Longitude, lat = ~Latitude, popup = content) %>%
  setView(-73.9949344, 40.7179112, zoom = 12) %>%
 addProviderTiles("OpenStreetMap.BlackAndWhite")

FireMap1

Question 2

Merged <- Merged %>%
  mutate(Prop_type = substr(PROPERTY_USE_DESC, 1, 1),
         Prop_type = substr(Prop_type, 1, 1),
         Prop_desc = ifelse(Prop_type == "2", "Educational", ifelse(Prop_type == "4", "Residential", ifelse(Prop_type == "5", "Commercial", ifelse(Prop_type == "6" | Prop_type == "7" | Prop_type == "8" | Prop_type == "9", "Industry/Infrastructure/Outdoor", "Other")))),
         Prop_desc = as.factor(Prop_desc))

pal = colorFactor("Set2", domain = Merged$Prop_desc) 
color_offsel1 = pal(Merged$Prop_desc)

content <- paste("Type of incident:",Merged$INCIDENT_TYPE_DESC,"<br/>",
                 "Action:",Merged$ACTION_TAKEN1_DESC,"<br/>",
                 "Detector Status:", Merged$DETECTOR_PRESENCE_DESC,"<br/>",
                 "Type of Property:", Merged$Prop_desc, "<br/r")

FireMap2a <- FireMap1 %>%
  addCircles(color = color_offsel1, popup = content, lng = ~Longitude, lat = ~Latitude) %>%
  addLegend(pal = pal, values = ~Merged$Prop_desc, title = "Type of Property")
  
FireMap2a
FireMap2a %>%
    addCircleMarkers(color = color_offsel1, popup = content, clusterOptions = markerClusterOptions())

Question 3

FireTruckIcon <- makeIcon(
  iconUrl = "http://www.clker.com/cliparts/P/P/1/P/t/C/fire-truck-md.png",
  iconWidth = 15, iconHeight = 15,
  iconAnchorX = 7.5, iconAnchorY = 8.5
)
FireMap2a %>%
  addCircleMarkers(color = color_offsel1, lng = ~Longitude, lat = ~Latitude, radius = Merged$UNITS_ONSCENE*0.5, fillOpacity=0.7, data = Merged, group = "Incidents") %>%
  addMarkers(data = Firehouses, group = "Firehouses", icon = FireTruckIcon) %>%
    addLayersControl(
    baseGroups = c("Incidents", 
                   "Firehouses"),
    overlayGroups = c("Firehouses", "Incidents"))

Question 4

mat <- distm(Merged[,c('Longitude','Latitude')], Firehouses[,c('Longitude','Latitude')], fun = distHaversine)

mat[is.na(mat)] <- 9999999

#Merged$FacilityName <- Firehouses$FacilityName[max.col(-mat)]
#Merged$Nearest_Distance <- apply(mat, 1, min)
Merged2 <- Merged %>%
  mutate(Nearest_Station = Firehouses$FacilityName[max.col(-mat)],
         Nearest_Distance = apply(mat, 1, min),
         Incident_time = mdy_hms(INCIDENT_DATE_TIME),
         Arrival_time = mdy_hms(ARRIVAL_DATE_TIME),
         Time_to_Arrive = (Arrival_time - Incident_time)/60 ,
         Residential = as.numeric(Prop_type == 4)) %>%
  filter(Nearest_Distance <= 2000 & Time_to_Arrive < 30)
## Warning in as.POSIXlt.POSIXct(x, tz): unknown timezone 'zone/tz/2018c.1.0/
## zoneinfo/America/New_York'
ggplot(data = Merged2, mapping = aes(x = Nearest_Distance, y = Time_to_Arrive, col = Prop_desc, group = Residential)) +
  geom_point(aes(x = Nearest_Distance, y = Time_to_Arrive, col = Prop_desc)) +
  ylab("Time to Arrival (Minutes)") +
  xlab("Distance to Nearest Station") +
  ggtitle("NYC Fire Response Times and Distance") +
  guides(col = guide_legend(title = "Property Type")) +
  geom_smooth(data = Merged2[Merged2$Prop_type == 4,], aes(y = Time_to_Arrive, x = Nearest_Distance), method = lm, se = FALSE, col = "Red") +
  geom_smooth(data = Merged2[Merged2$Prop_type != 4,], aes(y = Time_to_Arrive, x = Nearest_Distance), method = lm, se = FALSE, col = "Blue") +
  annotation_custom(grob = grid.text("Residential", x=0.9,  y=0.15, gp=gpar(col="Red", fontsize=9, fontface="bold"))) +
  annotation_custom(grob = grid.text("Non - Residential", x=0.85,  y=0.235, gp=gpar(col="Blue", fontsize=9, fontface="bold")))

## Don't know how to automatically pick scale for object of type difftime. Defaulting to continuous.

In order to better visualize the data, I eliminated some outliers both in distance (liklely due to miscoding the incident lat/long) and time (on incident took 60+ mins response time).

Overall we can see the expected pattern that incidents further from the closest fire station have a slightly longer wait between the onset of the incident and the arrival of the firs unit. There is another interesting pattern we can see from the plot above. In general the response time to residential units is lower than non-residential units. Moreover the the effect of distance on increased wait-time is quite a bit lower for residential units than non-residential units.

Merged3 <- Merged2 %>%
  group_by(ZIP_CODE) %>%
  mutate(avg_resp_time = mean(Time_to_Arrive)) %>%
  ungroup(ZIP_CODE) %>%
  select(ZIP_CODE, avg_resp_time) %>%
  rename(ZIPCODE = ZIP_CODE) %>%
  mutate(avg_resp_time = as.numeric(avg_resp_time))

Merged3 <- unique(Merged3)
nyc_census <- readOGR("Data/Polygon_NYC/.","ZIP_CODE_040114")
## OGR data source with driver: ESRI Shapefile 
## Source: "Data/Polygon_NYC/.", layer: "ZIP_CODE_040114"
## with 263 features
## It has 12 fields
nyc_census <- spTransform(nyc_census, 
                          CRS("+proj=longlat +datum=WGS84"))

nyc_census <- merge(nyc_census, Merged3, by = "ZIPCODE")
gg <- tm_shape(nyc_census) + layout +
tm_fill("avg_resp_time", title = "Average Response Time", n = 5, style = "quantile") +  # provides nicer intervals
tm_borders(alpha=.5)

gg

Based on the choropleth map above, we can see that there is a clear spatial heterogeneity when it comes to response times. Areas in Eastern Queens, Northern Bronx, and Southern Staten Island take have some of the longest wait times. Surprisingly a few areas of Manhattahn around midtown also have high wait times.